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Abstract 

Orbits  about  the  Martian  moons,  Phobos  and  Deimos,  previously 
found  to  be  stable  by  the  Poincar6  surface  of  section  technique  in  a 
restricted  three-body  system,  were  investigated  with  non-zero  eccentricity 
of  the  moons  added  to  the  dynamics.  Hamilton’s  canonical  equations  were 
derived  to  describe  the  dynamics  of  the  planar  four  body  system  (Mars, 
moon,  sun,  artificial  satellite).  Although  not  strictly  applicable  to  this  case, 
the  surface  of  section  technique  was  used  to  verify  whether  or  not  the  pre¬ 
viously  found  stable  orbits  were  realistic  when  a  real  world  perturbation 
was  added  to  the  system  dynamics.  The  surface  of  section  technique  in¬ 
volves  the  numerical  integration  of  several  orbits  with  the  same  initial 
value  of  the  Hamiltonian  based  on  matching  the  angular  rates  of  the  orbit 
and  the  moon  with  respect  to  Mars.  Apoapsis  and  periapsis  points  of  the 
orbits  were  plotted  in  the  two-dimensional  configuration  space.  Stable 
orbits  were  found  when  the  points  form  regions  which  are  bounded; 
chaotic  orbits  were  indicated  by  random  points. 
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REALISTIC  ORBITS 


ABOUT  THE  MARTIAN  MOONS 


T.  Introduction 

A  future  manned  mission  to  Mars  will  be  a  historic  achievement,  and 
the  possible  advancements  in  scientific  knowledge  of  the  solar  system  are 
unimaginable.  However,  when  the  overall  cost  of  such  a  venture  is  consid¬ 
ered,  space  exploration  is  also  astronomically  expensive. 

A  manned  mission  to  Mars  will  require  extensive  life  support 
equipment,  supplies,  and  repair  materials.  Basic  supplies  such  as  water, 
oxygen,  and  food  must  be  stored  in  the  spacecraft  with  the  astronauts  dur¬ 
ing  the  entire  length  of  the  mission.  Life  support  equipment  or  the  vehicle 
structure  may  need  repair  during  the  mission  which  necessitates  carrying 
along  tools  and  raw  materials.  These  items  must  also  be  stored  on  the 
spacecraft.  All  this  extra  weight  translates  directly  into  additional  cost  (in 
energy)  to  launch  the  spacecraft  from  earth. 

Due  to  the  earth’s  strong  gravitational  field,  any  extra  weight  is  very 
expensive  in  terms  of  the  amount  of  energy  (rocket  thrust)  needed  to 
escape  the  earth’s  gravity  and  enter  orbit.  If  earth’s  gravitational  field 
were  half  as  strong  there  would  be  a  significant  savings  in  energy  that 
could  be  used  elsewhere  in  the  mission  planning,  such  as  mission  duration 
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or  final  destination. 

Since  it  is  not  possible  to  change  the  earth’s  gravitational  field,  one 
logical  solution  is  to  find  a  small  moon  or  asteroid  that  can  act  as  a  way- 
station  or  refueling  stop.  This  way-station  will  have  a  smaller  gravitational 
field  than  the  earth’s  so  spacecraft  will  use  much  less  energy  to  land  and 
take  off.  While  on  the  way-station,  astronauts  could  mine  the  surface  for 
needed  materials,  minerals,  and  water  which  would  save  precious  energy 
and  extend  the  overall  duration  of  the  mission. 

For  a  manned  mission  to  Mars,  the  logical  choice  for  a  way-station  is 
one  of  its  small  moons.  Phobos  is  the  more  interesting  and  useful  choice 
due  to  its  larger  size  and  closer  proximity  to  Mars.  Both  moons  are 
believed  to  be  asteroids  that  were  captured  by  Mars’  gravitational  field. 
However,  their  extremely  small  size  and  close  proximity  to  Mars  means 
that  any  orbit  around  either  moon  will  be  greatly  influenced  by  the  gravity 
of  Mars  (1:18).  TTie  question  then  becomes  one  of  stability  of  the  orbit. 

When  the  Russians  launched  the  Phobos  2  scientific  probe  in  July 
1988,  it  was  the  first  Mars  mission  since  the  two  U.S.  Viking  probes  were 
launched  in  1975  (2:42).  The  data  from  the  Phobos  2  probe  answered 
many  questions  about  Mars  and  its  two  moons.  The  probe’s  unfortunate 
computer  failure  in  March  1989  did  not  diminish  the  success  of  its  57  days 
in  orbit  around  Mars  and  Phobos  (3:156). 

Updated  orbital  parameters  of  Phobos  will  allow  accurate  predictions 
of  the  moon’s  location.  This  new  data  reduces  the  uncertainty  of  the  moon’s 
position  from  500  kilometers  (km)  to  4  km  (2:46).  Since  the  diameter  of 
Phobos  is  only  27  km  on  its  longest  axis  (1:19),  a  500  km  error  in  predict- 
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ing  its  current  position  prevented  any  previous  rendezvous  or  landing 
attempt.  A  4  km  error  in  position  can  easily  be  handled  by  the  spacecraft’s 
computer  and  guidance  system. 

The  gravitational  field  of  Phobos  was  accurately  determined  before 
the  Phobos  2  probe  failed.  Escape  velocity  (speed  needed  to  escape  local 
gravity)  for  Phobos  is  8  miles  per  hour  (mph)  (1:18).  Earth’s  escape 
velocity  is  25,000  mph  and  requires  an  enormous  amount  of  highly-explo- 
sive  fuel  (1:17-18).  The  much  smaller  gravitational  field  on  Phobos  is  a 
major  reason  for  selecting  it  as  a  possible  future  way-station. 

Phobos  2  did  return  some  disappointing  data  about  the  amount  of 
water  on  the  surface  of  Phobos.  Although  many  scientists  want  to  believe 
Phobos  is  an  oasis  for  future  colonies  on  Mars,  the  Infrared  Spectrometer 
(ISM)  experiment  found  an  extremely  low  amount  of  hydration  in  the  first 
five  meters  of  the  surface  (3:158-159)  indicating  much  less  hydration  in  the 
regolith  (fractured  surface  layer)  than  on  most  known  asteroids  and  less 
than  what  has  been  found  on  the  surface  of  Mars  (3:159).  This  does  not 
eliminate  the  possibility  that  water  may  be  found  deeper  than  the  regolith. 
Even  if  Phobos  is  completely  dry,  it  is  still  an  excellent  choice  for  continu¬ 
ing  scientific  research  of  the  sun’s  gamma  rays  and  the  gravitational  and 
magnetic  fields  of  Mars. 
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n.  Aporoach/Methodology 


Previous  work  by  Jansson  (4)  proved  that  stable  orbits  about  the 
Martian  moons  are  possible,  but  the  use  of  several  simplifying  assumptions 
listed  below  may  make  the  orbits  unrealistic: 

Phobos  and  Deimos  have  circular  orbits  about  Mars. 

Mars  is  treated  as  a  point  mass. 

The  sun’s  gravity  field  is  ignored. 

Combining  these  assumptions  may  have  over-simplified  the  problem 
and  produced  unrealistic  solutions.  This  research  will  analyze  the  stable 
orbits  around  Phobos  and  Deimos  found  by  Jansson  to  see  if  the  orbits 
remain  stable  after  including  realistic  perturbation  terms  in  the  analysis. 

Stable  orbits  around  the  Martian  moons  will  be  investigated  using  the 
Poincar6  surface  of  section  stability  technique.  The  equations  of  motion 
will  be  derived  using  the  Hamiltonian  of  the  four  body  problem  (Mars, 
moon,  spacecraft,  sun).  Since  this  system  is  not  a  conservative  two  degree- 
of-freedom  Hamiltonian  system,  the  surface  of  section  technique  is  not 
strictly  applicable.  The  surface  of  section  stability  technique  numerically 
integrates  the  parameters  of  several  test  orbits  using  the  same  Hamiltonian 
value.  These  test  orbits  will  be  numerically  integrated  in  four-dimensional 
space  and  then  plotted  in  two-dimensional  space.  Each  dot  on  the  plot  will 
represent  a  closest  (periapsis)  or  farthest  (apoapsis)  point  in  the  orbit.  New 
dots  for  each  orbit  will  represent  a  small  perturbation  (change)  in  the 
orbit.  If  the  dots  produce  a  closed  circle  or  ellipse,  the  orbit  is  bounded 
and  therefore  stable.  Random  or  chaotic  patterns  will  indicate  unstable 
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orbits.  Plots  wiU  be  made  of  each  trajectory  to  see  if  the  orbit  remains 
stable  as  realistic  perturbations  are  added  to  the  system  dynamics. 
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m.  Problem  Dynamics 


The  coordinate  system  used  is  illustrated  in  Figure  1  showing  the 
moon  modeled  as  a  triaxial  ellipsoid  with  an  exaggerated  misalignment  of 
the  X  axis  due  to  the  eccentricity  of  the  moon. 

Jansson’s  research  assumed  the  moons’  orbits  around  Mars  were  cir¬ 
cular  (eccentricity  =  0).  Although  the  eccentricities  of  Phobos  and  Deimos 
are  small  (.015  and  .(XX)52  respectively)  (4:2),  this  investigation  will  not 
assume  circular  orbits  in  order  to  achieve  realistic  results.  Since  the  orbits 
will  not  be  circular,  the  constant  rotation  rates  (Q)  of  the  moons  about 
their  z  axis  will  not  be  equal  to  their  non-constant  revolution  rates  (“) 
about  Mars.  This  makes  the  dynamics  more  complicated  because  the 
moon’s  minimum  axis  of  inertia  (x)  will  not  remain  pointed  toward  Mars 
while  the  moon  rotates  about  its  maximum  axis  of  inertia  (z).  The  inter¬ 
mediate  axis  of  inertia  (y)  completes  the  orthogonal  set  of  axes. 

Also,  if  the  orbits  about  Phobos  continue  to  remain  stable,  the  J2 
term  of  the  geopotential  of  Mars  will  be  included.  This  investigation  will 
also  consider  the  sim’s  gravitational  attraction  modeled  as  a  point  mass  due 
to  its  distance  from  the  satellite. 

Equations  of  Motion 

The  equations  of  motion  of  the  satellite  can  be  derived  in  the  form  of 
the  Hamilton’s  canonical  equations.  The  Hamiltonian  is  defined  as 

H-IPiOi-L  (1) 
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Figure  1.  Coordinate  System 
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where  Pi  and  Qi  are  the  generalized  momenta  and  generalized  coordi¬ 
nates,  respectively.  L  is  defined  as  the  Lagrangian  and  is  found  by 


L  =  T- V 


(2) 


where  T  and  V  are  the  kinetic  energy  and  potential  energy,  respectively. 
The  first  step  in  deriving  the  equations  of  motion  is  to  define  T  and  V  of 
the  system  in  terms  of  the  generalized  coordinates  (  Oi 's  )  and  generalized 
velocities  (  Qi 's  ). 

Kinetic  Energy.  The  kinetic  energy  of  the  satellite  as  it  orbits  about 
the  Martian  moon  is  given  by 


(3) 


where  m  is  the  mass  of  the  satellite  and  v  is  the  inertial  velocity  of  the 
satellite.  Based  on  very  small  moon  masses,  the  the  inertial  reference 
frame  is  set  at  the  mass  center  of  Mars. 

Now 


V  =  d  (4) 

where  a  dot  over  a  vector  represents  the  time  rate  of  change  of  the  vector 
with  respect  to  the  inertial  frame  (inertial  derivative). 
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Therefore, 


and  from  Figure  1 


V  ^  =  d  •  d 


d  =  D  +  R 


leading  directly  to 


d  =  5  +  R. 


(5) 


(6) 


(7) 


Here,  R  may  be  written  as 


(8) 


where  the  superscripts  i  and  r  indicate  derivatives  taken  in  the  inertial  and 
rotating  (moon-fixed)  frames,  respectively,  with  the  moon  rotating  at  a 
fixed  rate. 

From  Figure  1, 


and 


ri  -•  ^ 

(o  =  Q  =  Q  az 


(9) 


R  —  X  Y  ^  Z  Ax 


(10) 
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which  yields 


[(nSJxl 


X  3x  Y  Z  3.z)j 


(11) 


and  simplifies  to 


R  =  (x  -  Qy}  ax  +  (y  +  nx)  ay  +  Z  az 


(12) 


Since  we  are  not  assuming  the  moon  is  in  a  circular  orbit  about  Mars,  the 
D  term  becomes  complicated  because  the  moon’s  X  axis  does  not  remain 
fixed  on  the  mass  center  of  Mars  as  the  moon  revolves  around  Mars. 
Taking  this  into  account,  the  magnitude  of  D  will  not  remain  constant  and 
will  be  calculated  by  (5:76) 


D-  a(l  -e^) 

1  +  e  CCS  V  ( 1 : 

where  the  true  anomaly,  v,  is  defined  as  the  angle  from  the  moon’s  closest 
approach  (perimartem)  to  the  current  position  vector  D.  The  mean 
anomaly  (M)  is  the  mean  motion  (n)  of  the  moon  multiplied  by  the  time 
since  perimartem  passage.  The  difference  between  the  moon’s  x  axis  and 
D  from  the  mass  center  of  Mars  is  the  difference  between  v  and  M.  This 
is  sometimes  called  the  “equation  of  center”  (5:77)  and  is  shown  in 
Figure  1  as  v  -  M.  The  letters  a  and  e  refer  to  the  moon’s  orbital  semi¬ 
major  axis  and  eccentricity,  respectively. 

Approximating  equation  (13)  by  a  periodic  expansion  in  e  and  dis- 
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carding  fourth  order  terms  and  higher  (5:76),  because  their  contribution  is 
negligible  (<  .5  km),  yields 


D»a(l  +  —  e^-e  cos  M  -  —  e  ^  cos  2M  +  —  e  ^  cos  M 
^2  2  8 


-  —  e  ^  cos  3M  H — ) 


Similarly,  the  “equation  of  center”,  defined  as  v  -  M,  can  be  approximated 
by  (5:77) 


(v  -  m)  »  2c  -  —  sin  M  +  —  e^  sin  2M  +  e^  sin  3M  h — 

\  4  /  4  12 


The  position  vector  D  can  now  be  written  as 


D  =  D  cos(v  -  m)  ax  -  D  sin(v  -  m)  ay 


As  in  equation  (8),  D  can  be  written  as 

‘|5  =  'AD.t"xD) 


The  first  term  on  the  right  hand  side  of  equation  (17)  is  written  as 


D  =  cos  (v  -  m)]  ax  -  [d  sin  (v  -  m)]  aj 


(18) 


and  the  product  rule  yields 


D  =  [d'  cos  (v  -  m)  -  D  (v  -  Mf  sin  (v  -  m)]  ax 
-  [d*  sin  (v  -  m)  +  D  (v  -  Mf  cos  (v  -  m)] 


where 


=  aQ  (e  sin  M  +  e^  sin  2M  -  —  e^  sin  M  +  —  sin  3M 

8  8 


(v  -  m)  =  O  f2e  -  -  e^l  cos  M  +  —  cos  2M  +  —  cos  3M 
4/2  4 


The  second  term  on  the  right  hand  side  of  equation  (17)  is  not  as 
complicated  and  is  written  as 

CO  X  D  =  Q  az  X  (d  cos  (v  -  m)  ax  -  D  sin  (v  -  m)  ay)  i 

and  yields 


(0  X  D  =  nD(  sin  (v  -  m)  ax  +  cos  (v  -  m)  ay] 


Combining  equations  (19)  and  (23)  and  then  simplifying  results  in 


D  = 


D  cos  (v  -  m)  +  D  sin  (v  -  m)  |q  -  (v  -  M)'j 
D  cos  (v  -  m)  |q  -  (v  -  m)')  -  D  sin  (v  -  m) 


ax 

a 


y. 


(24) 


Combining  equations  (12)  and  (24)  into  equation  (7)  gives 

d  =  X  -  QY  +  D  cos  (v  -  m)  +  D  sin  (v  -  m)  jo  -  (v  -  m)' ) 
+  Y  +  nx  +  D  cos  (v  -  m)  (q  -  (v  -  m)' )  -  D'sin  (v  -  m) 

+  2  32 


ax 


(25) 


Since  the  dynamics  are  assumed  to  remain  in  the  same  plane  of  motion,  the 
Z  term  will  be  set  to  zero  so  the  surface  of  section  technique  can  be  used. 

Finally,  combining  equations  (25),  (5),  and  (3)  gives  the  total  kinetic 
energy  of  the  satellite  as 

Tsat  =  j  m„,  [x  -  QY  +  D  cos  (v  -  m)  +  D  sin  (v  -  m)  (q  -  (v  -  m)  )]  ^ 

+  -  m,.,  [y  +  nx  +  D  cos  (v  -  m)  (n  -  (v  -  m)  )  -  D'sin  (v  -  m)]  ^ 

2  (26) 


Potential  Energy.  The  potential  energy  of  the  satellite  can  be  divided 
into  three  parts 


Vsat  —  V moon  V Mars  ^ sun 


(27) 


Potential  Energy  Due  to  the  Moon.  The  potential  energy  of 
the  satellite  due  to  the  gravitational  attraction  of  the  moon  is  not  simple  due 
to  the  non-spherical  shape  of  the  moon.  'Fhe  moon  can  be  modeled  as  a  tri- 
axial  ellipsoid  as  shown  in  Figure  2.  The  gravitational  potential  for  an 
arbitrarily  shaped  body  is  given  by  Meirovitch  (6:430-436)  starting  with 
the  integral  form  of  the  inverse  square  force  field  formula 


(28) 

where  G  is  the  universal  gravitational  constant  and  from  Figure  2 


V  moon  — Gmsat 


I 


dM 


moon 


r  =  |;|=  R  +  p| 

(29) 

p  =  a  ax  +  b  ay  c  az 

(30) 

p^  =  a^ 

(31) 

r2  =  x2  +  y2  +  z2 

(32) 

R  •  p  =  aX  +  bY  +  cZ 

(33) 

Now  equation  (29)  can  be  moved  to  the  numerator  and  written  as 


r‘  = 


R-P 


=[(R-pr] 


1/2 


r-l 


=  (r2-2R  •  P  +p")' 


(34) 

(35) 


1  4 


^2 


ax 


Figure  2.  Gravitational  Potential  Due  to  Moon 
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-  1/2 


r 


-1 


2R  •  p 
R2 


r-1  =  R-i 


2R  •  3 

1 - 

R2 


(36) 


(37) 


Using  a  binomial  expansion  and  ignoring  terms  higher  than  second  order 
because  their  contribution  is  negligible  (<  .5  km),  equation  (37)  can  be 
written  in  the  form  of  a  power  series  as 


r 


-1 


1  + 


+ 


^2^ 
8\  r2 


r-l  =  R-l  +  R-3(r  .  p)  -  iR-3p2  +  |r-5|r  .  p)^+ •  •  • 


(38) 

(39) 


Combining  equations  (28),  (31)-(33)  and  (39)  yields 


moon 


■ 


*  dMmoon  -  (aX  +  bY  +  cZ)  dM^oon 

R3 


(a^  +  +  c^)  dMmoon 

2R3  J 


3Gm 


2R5 


(aX  +  bY  +  cZ)  2  dMmoon 


1  6 


(40) 


Integrating  the  first  term  of  equation  (40)  is  trivial  and  the  second 
term  vanishes  because  the  origin  is  taken  to  be  at  the  center  of  mass  leaving 
the  last  two  terms  which  can  be  combined  as 


[3R-2(aX  +  bY  +  cZp-  (a2+b2+c2)]  dMr 

2R3 


-Gmsat 

2R3 


{3R-2x2  -  l)|  a2dMmoon-^^3R-2Y2  -  l)|  b^dMmoon 
-43R-2Z2  -  l)|  c^dMmoon 

+  6R-2|xY  j  abdMmoon  +  Xz|  acdMmoon  +  Yz|  bcdMp 


Defining  the  mass  moments  and  mass  products  of  inertia  from  Likens 
(7:525),  it  follows  from  Meirovitch  (6:435)  that 


j  a^  dMmoon  —  ~  (Izz  lyy  "  Ixx) 


b^  dMmoon  —  ~  (Izz  +  Ixx  "  lyy) 
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I 


c2  dM 


moon 


=  ^{I 
2 


XX  lyy  "  ^zz) 


(45) 


I 

I 

I 


ab  dMmoon  —  Ixy 


ac  dMmoon  —  Ixz 


be  dMmoon  —  lyz 


(46) 


(47) 


(48) 


using  the  principal  axes  frame,  we  have 

Ixy  =  Ixz  =  lyz  =  0. 


(49) 


This  simplifies  equation  (42)  allowing  the  total  potential  energy  of  the 
satellite  due  to  the  moon  to  be  written  as 


moon 


GmsatMmoon  Gnisat 
R  4r3 


{Ixx  +  lyy  +  Izz) 


3Cjmsat 

4R5 


[X^lyy  +  I 


zz 


SCjiTisat 

4R5 


[y^Ixx  Izz  *  lyy)] 


4R5 


yy 


(50) 


Potential  Energy  Due  to  Mars.  The  potential  energy  of  the 
satellite  due  to  the  gravitational  attraction  of  Mars  is  the  Newtonian  term 
plus  the  J2  term  of  the  geopotential  simplified  due  to  the  planar  restric¬ 
tions. 


^Mars  — 


GmsatMMars 

d 


1  + 


^  d 


Now  d  is  defined  as 


(51) 


d  =|d|  =|D-t-R 


(52) 


and  D  and  R  are  defined  by  equations  (10)  and  (16)  respectively  and 
yield 


d  =  |[x  +  D  cos  (v  -  m)]  ax  +  [y  -  D  sin  (v  -  m)]  ay 


d  =  V[X  +  D  cos  (v  -  m)]  ^  +  [y  -  D  sin  (v  -  m)]  ^ 


(53) 

(54) 


Potential  Energy  Due  to  the  Sun.  The  potential  energy  of  the 
satellite  due  to  the  gravitational  attraction  of  the  sun  can  be  treated  as  a 
third  body  perturbation  (8:89-95).  The  disturbing  function  for  general 
third  body  perturbations  can  be  determined  from  Figure  3  and  is  written  as 


Vj  —  Gmsat  Mj 


L.^ 

\P2j 


PJ 


(55) 


1  9 


Figure  3.  Coordinate  System  for  General  Third  Body 
Disturbing  Function 


20 


For  the  sun-Mars-satellite  system  it  may  be  written  as 


V  sun  —  Gmsat^sun 


1 


iP2sun 


Psun  ■  d 
3 

Psun  / 


(56) 


where  P2sun  is  the  vector  from  the  satellite  to  the  perturbing  mass  (the 
sun)  and  Psun  locates  the  third  mass  relative  to  the  primary  body  (Mars). 
The  first  term  on  the  right  hand  side  of  equation  (56)  is  called  the  direct 
term  and  represents  the  direct  gravitational  attraction  of  the  third  mass  on 
the  satellite.  The  second  term  on  the  right  hand  side  of  equation  (56)  is  the 
indirect  term  and  represents  perturbations  of  the  third  mass  on  the  primary 
body. 

Now  the  indirect  term  can  be  written  as 


V  indirect  —  G^saiMsun 


- 

Psun  *  d  1 

—  GMjun 

1  ^  \ 

3  I 

2 

Psun  / 

iPsun/ 

I  cos  S 


(57) 


and  the  direct  term  becomes 


Vdirect  =  -  (l  .  2acos  S  +  a^)' 

P2sun  Psun 


(58) 


letting 


Psun 


2  1 


(58.1) 


If  the  ratio  ot  is  quite  small  (as  it  is  in  this  case),  a  direct  expansion 
in  powers  of  ct  will  be  rapidly  convergent.  This  approach  is  called  the 
Lunar  Theory  (8:91)  and  allows  a  binomial  expansion  of  equation  (58) 
which  yields 


Vdirect  —  “ 


sati^^sun 


Psun 


1  +  a  cos  S  +  a^|-^  +  |(cos  S)^j  +  • 


(59) 


The  first  term  in  this  expansion  does  not  contain  the  coordinates  of 
the  satellite  and  therefore  can  be  ignored.  The  second  term  exactly  cancels 
the  ^indirect  term  leaving 


v,„„= 

2psun 

where  (cos  S)  can  be  written  as 


(60) 


cos  S  = 


Piun  •  R 
Piun  |R 


(60.1) 


with  R  defined  by  equation  (10)  and  Psun  will  be  approximated  by 


Psun  ”  p<iinCOS  M  Ss  "  p«i»iSin  M  Sy 


(60.2) 
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Finally,  equation  (27)  can  now  be  written  as 


Vsat  — 


GnisatMmoon  Gmsat  ,  .  t  \ 

+ - r—  (Ixx  *yy  ^zzf 


R 

3Cjinsat 
4R5 
3Cjinsat 
4R5 
CjmsatMMarsI 


4R3 

[X2(Iyy+I^^-l^j 

(Ixx  Izz  ■  lyy)] 

J  ^J2/Rm«.P 


d  I  2 

Gnisat^sund^ 


^  3 
2Psun 


{3(cos  Sf-l) 


(61) 


Lagrangian.  The  Lagrangian  is  now  needed  and  is  formed  by  com¬ 
bining  equations  (2),  (26),  and  (61).  Dividing  out  the  mass  of  the  satellite 
allows  the  Lagrangian  to  be  written  in  a  per  unit  mass  basis  as 


=  i[x-i:2Y-fD'  cos(v-M)  +  Dsin(v-M)(Q-(v-M)| 
+  j  [ Y  +  QX  -  D'  sin(v  -  m)  +  D  cos(v-  m)(q  -  (v  -  m)' 

3G 


4R5 

3G 


(lyy  +  Izz  ■  Ixx)] 

■  +  Izz  '  lyy)] 

GMMara  ,  J2/RM  arsP 

dH  d  I 
+  5Mffi!!£(3(co5S)2-l) 


^Psun 


(62) 
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Generalized  Momenta.  The  generalized  momenta  are  derived  as 


Px  =  — =  X-nY  +  D'cos(v-M)  +  Dsin(v-M)(n-(v-M)') 

(63) 

Py  =  —  =  Y  +  QX  -  D'  sin(v-  m)  +  D  cos(v  -  m)(q  -  (v  -  m)1 

(64) 


When  compared  to  equation  (25),  it  is  useful  to  note  that  the  generalized 
momenta  equations  are  the  inertial  velocity  components  resolved  on  the 
moon  axes  and  will  be  used  later  to  determine  initial  momenta  values. 


Generalized  Velocities.  Now  the  generalized  velocities  (Oi)  may  be 
eliminated  from  the  Lagrangian.  Equations  (63)  and  (64)  are  rearranged 
to  solve  for  X  and  Y 


X  =  Px  n Y  -  D'  cos  (v  -  m)  -  D  sin  (v  -  m)  (q  -  (v  -  m)') 


(65) 


Y  =  Py  -  QX  +  D'  sin(v - m)-  D cos(v  -  M)(£i  -  (v  -  mJ) 


(66) 


which  are  then  substituted  back  into  the  Lagrangian  to  give 


3G 


[x^  (lyy  +  Izz  -  Izx)]  (Ixx  +  Izz  “  lyy)] 


GMMars 


1  + 


j  z/Rmi 


fl 


GMsund^ 

2Psun 


(3  (cos  Sf  -  \) 


(67) 
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The  Lagrangian  is  now  in  the  correct  form  to  use  equation  (1). 


1.  From  equation  (1)  the  Hamiltonian  can  be  written  as 


=  Px(Px  +nY-D'cos(v-M)-Dsin(v-M)(Q-(v-M 
t-  Py(py-  iiX  +  D'sin(v-M)-Dcos(v-M)(n-(v-M)) 


1  [p2  p2]  GMm 

2[Px  +  PyJ-  r 


moon 


(Ixx  +  lyy  fzz) 


^  [X2  (lyy  +  -  I^x)]  -  ^  [  Y2  (Ixx  +  Izz  -  lyy)] 

4R^  4R^ 


GMMars 

d 


1  +  1  .  (3  SF  -  l) 

2\  d  } \  ^2 

■^Ksun 


or  with  some  simplification  as 


H  =  Px(nY  -  D'  cos  (v  -  m)  -  D  sin  (v  -  m)  (q  -  (v  -  m)  )) 

+  Py  (  -  QX  +  D'  sin(v  - m)-  D cos(v  -  M)(ft - (v  -  m)D 

-*■  |[Px Py]  -  ^yy 

-  [x^  (lyy  +  Izz  -  Ixx)]  "  [  Y^  (Ixx  Izz  '  lyy)] 


GMMaTiJ  ,  J2  Rm« 
d  2d 


GMsund^ 

2Psun 


(3  (cos  Sf  -  1) 
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Hamilton's  Equations.  The  equations  of  motion  can  now  be  written 
in  canonical  form  which  will  define  the  orbital  characteristics  of  the  sys¬ 
tem. 

X  =  ^  =  Px  +QY  -  D'cos(v-M)-Dsin(v-M)(Q-(v-M)') 

dPx  (70) 

Y  =  =  Py  -  nx  +  D'  sin (v-  m)-  D  cos (v -  m)(q  - (v  -  m)') 

dPv  nn 


f.  ^  GMmoonX  .  3GX  ,  , 

=  'dX  “  - ^5 - +  ^(3Iyy  +  3Izz-Ixx) 


dx 

15GX 


[X2  (lyy  -K 

Izz  -  Ixx)  +  Y^  (Ixx  +  Izz  ■  lyy)] 


GMMars(X-»-Dcos(v-M))  ^  ^  3J2R^„, 

[  2d2 

GMsurjx -t- Dcos(v-M))  I  3A  A 

p3  (X2  +  y2)  ■ 

Psun 

GMsun3d^  (x(cos  M)f^  -  Ysin  Mcos  m)  XA 

pL  L  '(X^  +  Yi)^ 


dH 

dY  “  ■ 
15GY 


_  r»  /-Y  GMmoonY  3GY  ^  ^  ^ 

=  -  Pxi2 - ; - + - r  {3Ixx  +  3Izz  -  lyy) 

R3  4R5 


[X2  (lyy  +  Izz  ■  Ixx)  (Ixx  +  Izz  ■  lyy)] 


GMMarjY-Dsin(v-Mtf  ^  3J2R^„, 
d^  [  2d^ 

GMsun(  Y  -  Dsin  (v-m))  /  3A  .  \ 

Psun 

GMsun3d^  (yisin  M)^  -  Xsin  Mcos  m)  YA 
dL.  (X2  +  Y2)  (x2  +  Y2)2 
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letting  A  be  defined  as 


A  =  (Xcos  M  -  Ysin  Mf 


(73.1) 


X-Axis  Symmetry 

Jansson  found  that  the  equations  of  motion  may  be  transformed  using 
the  following  relationships 

X  =>  X 

y  =>-y 

t  =>  -t 

without  effecting  the  equations  of  motion  in  the  rotating  coordinate  system 
of  the  restricted  three  body  problem.  The  result  is  a  second  orbit  which  is 
a  reflection  of  the  first  about  the  x  axis  and  traversing  in  the  opposite 
direction  from  the  first  orbit.  This  effectively  doubles  the  number  of  data 
points  obtained  for  a  given  trajectory  and  is  used  to  duplicate  Jansson ’s  sur¬ 
face  of  section  plots.  However,  it  will  not  be  used  after  perturbations  are 
added  to  the  system  since  perturbations  destroy  the  symmetry. 
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IV.  Surface  of  Section  Technique 


In  this  problem  there  are  four  initial  conditions  defined,  and  if  an 
integral  of  motion  exists,  solutions  of  the  system  may  be  represented  in  a 
three-dimensional  subspace  for  a  given  value  of  the  constant  integral  of 
motion  (9:127).  This  case  is  not  a  restricted  three-body  problem  because 
the  eccentricity  of  the  moon  is  included.  Therefore  the  surface  of  section 
technique  is  not  strictly  applicable  but  is  used  in  an  attempt  to  verify  if  the 
periodic  orbits  found  by  Jansson  are  realistic  and  remain  stable. 

Poincare’s  surface  of  section  technique  is  an  analytical  technique 
which  is  widely  used  in  numerical  exploration.  The  restricted  three-body 
problem  has  a  four  dimensional  phase  space  and  one  known  integral  of 
motion.  By  considering  all  orbits  with  a  fixed  value  of  the  Hamiltonian  the 
resulting  group  of  orbits  must  fill  a  three  dimensional  subspace  of  the  full 
four  dimensional  phase  space.  If  a  second  integral  of  motion  exists,  it 
would  further  restrict  the  motion  of  the  orbits  and  constrain  each  orbit  to 
lie  on  a  two  dimensional  surface  imbedded  in  the  three  dimensional  sub¬ 
space.  If  a  plane  (called  the  surface  of  section)  is  passed  through  this  three 
dimensional  subspace  in  any  arbitrary  way,  the  internal  structure  of  the 
three  dimensional  subspace  can  be  examined  in  two  dimensions.  When  a 
trajectory  on  the  three  dimensional  manifold  passes  through  the  surface  of 
section  plane  a  point  of  crossing  can  be  plotted.  As  more  points  are  plot¬ 
ted,  stable  periodic  orbits  will  be  shown  as  a  series  of  closed  curves.  If  a 
second  integral  of  motion  does  not  exist,  the  trajectories  can  intertwine 
randomly  in  the  three  dimensional  subspace  eventually  filling  it  up  and  the 
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plotted  points  on  the  surface  of  section  will  be  chaotic. (8: 128) 

This  technique  is  used  to  examine  periodic  orbits  around  the  Martian 
moons  where  the  surface  of  section  is  defined  by  the  following  condition 

R  •  V  =  0  (74) 

or  for  the  planar  case  under  consideration 

XX  +  YY  =  0  (75) 

This  condition  occurs  at  the  satellite’s  apoapsis  and  periapsis  points.  The 
points  where  the  satellite  is  farthest  from  and  closest  to  the  moon.  These 
points  are  plotted  in  the  two  dimensional  configuration  space. 

Other  surface  of  section  criteria  are  possible  but  this  one  immedi¬ 
ately  indicates  if  the  trajectory  will  collide  with  the  moon.  If  any  one  of 
the  surface  of  section  points  has  a  magnitude  less  than  that  of  the  moon,  it 
will  obviously  collide  with  the  moon  before  completing  the  orbit. 

Stability  of  Orbits 

For  a  given  Hamiltonian  value,  stable  orbits  will  be  found  if  the  sur¬ 
face  of  section  points  form  closed  contours.  A  periodic  orbit  will  pass 
through  the  centers  of  the  closed  curves. 

If  the  closed  curves  intersect  one  another  (unstable  equilibrium  point 
or  saddle  point),  unstable  orbits  may  be  located  which  pass  through  the 
intersection  points. 
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V.  Solution  Method 


A  FORTRAN  computer  program  was  written  to  integrate  the  equa¬ 
tions  of  motion  using  the  Haming  fourth  order  integration  algorithm  to 
produce  satellite  trajectories.  Initial  satellite  X  and  Y  coordinates  are  input 
along  with  the  Hamiltonian.  The  output  is  a  data  file  with  the  surface  of 
section  X  and  Y  coordinates  that  can  be  plotted  by  a  two-dimensional  plot¬ 
ting  program. 


Initial  Conditions 

Initial  conditions  for  the  state  vector  (X,  Y,  Px,  Py)  must  be  deter¬ 
mined.  The  initial  values  of  X,  Y,  and  the  Hamiltonian  are  chosen  so  that 
the  surface  of  section  conditions  (equations  74  and  75)  are  met.  Given  this, 
the  initial  values  of  Px  and  Py  may  be  calculated  by  first  obtaining  an 
expression  for  Y. 

Substituting  equations  (63)  and  (64)  into  the  Hamiltonian  (69)  and 
combining  all  the  potential  energy  terms  into  the  quantity  V  results  in  the 
following 


H  =  |[x^+Y^  -  Q^X2+Y2)  -  D'^  -  (d(q  -  (v-M)')f  ] 

+  noia  -  {v-M)tYsin  (v-M)  -  Xcos  (v-M)] 

+  QD^Xsin  (v-m)  +  Ycos  (v-m)]  +  V  (76) 


Solving  equation  (75)  for  X  gives 


(77) 
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Substituting  equation  (77)  into  equation  (76)  and  rearranging  yields 


(X2  +  Y2)  =  2X2  H  -  V  +  y2)^.  +i(D(n  -  (v-M)')f 

-  2x200(0  -  (v-Mltvsin  (v-M)  -  Xcos  (v-m|] 

-  2X2QDf  Xsin  (v-m)  +  Ycos  (v-M)] 


(78) 


Solving  equation  (78)  for  Y  results  in  two  values  which  Jansson  found  to 
yield  the  same  trajectory  but  traveling  in  opposite  directions  (prograde  and 
retrograde)  for  the  restricted  three-body  problem  with  zero  eccentricity. 

X  may  be  calculated  by  equation  (75)  and  then  the  required  initial 
values  for  Px  and  Py  may  be  determined  from  equations  (63)  and  (64). 

Trajectory  Integration 

Haming’s  Ordinary  Differential  Equation  Integrator  is  used  to 
numerically  integrate  the  equations  of  motion.  Haming’s  code  is  a  fourth 
order  predictor-corrector  algorithm  which  is  numerically  very  stable  and 
provides  high  accuracy.  It  extrapolates  the  last  four  values  of  the  state 
vector,  predicts  the  next  value,  evaluates  the  equations  of  motion  at  the 
predicted  value,  and  then  corrects  the  value  using  a  higher  order  poly¬ 
nomial  (8:120). 

Jansson  was  able  to  use  an  integration  step  size  of  200  seconds  for 
most  trajectories  but,  after  adding  in  the  eccentricity  to  the  dynamics,  an 
integration  step  size  of  50-100  seconds  was  required  to  prevent  numerical 


integration  problems.  Even  with  these  large  integration  step  sizes, 
extremely  accurate  surface  of  section  points  were  determined  using  an 
interpolating  polynomial. 

Surface  of  Section  Points 

At  each  integration  step,  the  current  value  of  R  •  V  is  compared  to 
the  last  and  if  the  sign  has  changed,  the  surface  of  section  condition  stated 
in  equation  (74)  has  been  met. 

To  obtain  accurate  values  for  the  surface  of  section  points,  a  Newton 
Interpolation  Polynomial  (10:112)  is  used  with  the  last  four  values  inte¬ 
grated.  Since  R  •  V  is  a  function  of  X  and  Y,  inverse  interpolation 
(10: 1 19)  is  used  to  determine  the  X  and  Y  coordinates  for  which  R  V  =  0, 
Remembering  that  there  is  trajectory  symmetry  about  the  X  axis  (for 
Jansson’s  problem)  gives  a  reflected  surface  of  section  point  (X,  -Y)  for 
every  point  (X,  Y)  obtained  from  the  integration.  There  will  be  no 
reflected  surface  of  section  points  used  once  eccentricity  is  added. 

Orbit  Checks 

Checks  are  made  throughout  the  integration  to  ensure  the  trajectory 
has  not  escaped  the  moon’s  gravitational  influence  or  passed  too  closely  to 
the  mass  center  of  the  moon  which  may  create  a  singularity  in  the 
integration.  If  either  of  these  occur,  the  integration  is  halted. 

Verification  of  Dynamics  Model  and  Computer  Code 

Several  checks  were  made  to  verily  that  the  equations  of  motion  and 
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the  FORTRAN  computer  code  were  correct.  One  check  was  to  observe  the 
Hamiltonian.  The  Hamiltonian  is  a  constant  integral  of  motion  and  should 
remain  invariant  (to  several  decimal  places)  over  the  entire  trajectory  and 
was  checked  periodically  throughout  the  integration  to  insure  that  it  was 
conserved  while  eccentricity  was  kept  equal  to  zero. 

Another  check  was  to  reduce  the  computer  code  dynamics  to  the 
simple  two-body  problem.  To  do  this,  the  gravitational  attraction  of  Mars 
and  the  sun,  the  system  rotation,  the  moon’s  eccentricity  and  mass  moments 
of  inertia,  and  the  distance  from  the  moon  to  the  origin  of  the  coordinate 
system  were  all  set  to  zero.  The  satellite  was  then  given  circular  velocity 
based  on  an  initial  position  relative  to  the  moon.  As  expected,  the 
integrated  trajectory  was  circular,  returned  to  the  original  starting 
position,  and  maintained  a  constant  Hamiltonian.  When  the  moon’s  mass 
moments  of  inertia  were  included,  the  orbit  was  no  longer  circular  but  the 
Hamiltonian  did  remain  constant  to  eight  decimal  places. 

The  next  check  involved  setting  the  gravitational  attraction  of  the 
moon  and  the  sun,  the  ^2  term  of  the  geopotential  of  Mars,  and  the  initial 
velocity  of  the  satellite  to  zero.  The  satellite  position  was  chosen  to  lead 
the  moon  by  a  small  distance.  The  integrated  trajectory  showed  the  satel¬ 
lite  remained  relatively  stationary  in  the  rotating  coordinate  system,  indi¬ 
cating  is  was  in  the  expected  circular  orbit  about  Mars. 

Next,  previous  research  output  should  be  reproducible  given  the 
same  equations  of  motion  and  initial  conditions.  To  duplicate  Jansson’s 
wofk,  the  gravitational  attraction  of  the  sun,  the  ^2  term  of  the  geopotential 
of  Mars,  and  the  eccentricity  of  the  moon  about  Mars  must  all  be  set  to 
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zero.  A  surface  of  section  plot  can  then  be  duplicated  from  his  initial 
conditions. 

As  expected,  the  surface  of  section  plots  from  Jansson  were  repro¬ 
duced  for  both  Phobos  and  Deimos.  Typical  surface  of  section  plots  for 
Phobos  and  Deimos  are  compared  to  Jansson ’s  work  with  the  same  initial 
conditions  and  are  shown  in  Figures  4-7. 

Initial  Positions 

For  all  integration  mns  the  initial  position  of  the  satellite  was  set  at 
(-X,  0)  with  values  of  X  matching  Jansson ’s  periodic  orbits.  The  initial 
position  chosen  results  in  both  the  initial  velocity  and  the  initial  momenta  in 
the  X  direction  being  equal  to  zero. 
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Figure  5.  Janssen’s  Surface  of  Section  Plot  for  Phobos, 
H  =  -6.8528  (4:32) 
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Figure  6.  Recreation  of  Dcimos  Surface  of  Section  Plot, 

H  =  -2.738592 
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Figure  7.  Jansson’s  Surface  of  Section  Plot  for  Deimos, 
H  =  -2.738592  (4:52) 
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VI.  Results 


With  the  computer  model  and  dynamics  verified  the  non-circular 
eccentricity  of  the  moon  was  introduced  to  the  system.  Although  the 
eccentricities  of  the  moons  are  small,  they  cannot  be  ignored  for  any  seri¬ 
ous  orbital  predictions. 

As  long  as  the  eccentricity  was  zero,  the  Hamiltonian  could  be  used 
to  calculate  the  initial  velocity  and  therefore  the  initial  momenta  values  as 
discussed  previously.  But  the  Hamiltonian  will  not  remain  constant  once 
the  eccentricity  has  been  introduced  which  necessitates  a  new  method  of 
determining  the  initial  momenta  values.  The  first  approach  tried  used  the 
initial  velocities  as  calculated  from  Jansson’s  quadratic  equation  (4:24). 
This  value  is  identical  to  equation  (78)  when  the  initial  conditions  are 
applied.  Using  this  initial  velocity  proved  inappropriate  for  this  problem 
and  resulted  in  the  immediate  escape  of  all  trajectories  from  both  Phobos 
and  Ueimos.  Figure  8  is  a  plot  of  a  typical  trajectory  using  this  initial 
velocity  and  clearly  shows  an  eccentric  orbit  about  Mars  which  is  not 
tracking  the  revolution  rate  of  the  moon  about  Mars.  This  indicates  that  a 
much  more  refined  initial  condition  is  required. 

The  momenta  can  be  calculated  with  a  forced  value  of  the  orbital 
period  which  will  match  the  rates  of  the  orbit  to  that  of  the  moon  with 
respect  to  Mars.  In  the  inertial  reference  frame. 


d  =  Pxax  +  Pyay 


(79) 
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as  defined  by  equations  (25)  and  (63). 

Using  equation  (5)  along  with  the  previously  stated  initial  conditions 
leads  directly  to 


v2  =  pj 


(80) 


Now  using  the  period  of  the  moon’s  orbit  about  Mars  (1 1:33) 


TP  = 


2nV? 


(81) 


solving  for  the  semi-major  axis  gives 


a  = 


VITPPGMm, 


ars 


(2nf 


(82) 


but  a  can  also  be  written  in  terms  of  the  specific  mechanical  energy  (11 :28) 


£  =  - 


OMmot 

2a 


(83) 


which  then  gives 


Py  =  ± 


(85) 


Since  the  initial  position  is  (-X,  0)  and  the  orbits  being  investigated  are  all 
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retrograde,  the  only  possible  root  for  the  momenta  is  the  positive  one. 

The  momenta  for  each  periodic  orbit  found  by  Jansson  (4:33,50)  are 
calculated  in  the  manner  described  above  and  are  then  fixed  in  the 
FORTRAN  program.  Tables  1  and  2  show  the  results  for  Phobos  and 
Deimos  respectively.  The  FORTRAN  program  then  uses  the  fixed  values 
of  the  momenta  to  initialize  the  state  vector. 

The  periodic  orbits  of  the  most  practical  interest  for  use  as  parking 
orbits  are  the  inner  orbits,  that  is  the  orbits  which  stay  the  closest  to  the 
moon.  The  inner  orbits  did  not  handle  the  eccentricity  perturbation  very 
well  and  either  collided  with  the  surface  of  the  moon  or  escaped  into  an 
eccentric  orbit  about  Mars.  This  leads  us  to  the  conclusion  that  our 
assumption  that  the  orbits  are  reaUy  in  orbit  about  Mars  does  not  hold  close 
in  to  the  moon.  Results  for  each  moon  will  be  discussed  separately. 

OlQbQS 

With  only  the  eccentricity  of  Phobos  added  to  the  system  dynamics, 
the  five  closest  periodic  orbits  all  collided  with  Phobos  ranging  fmm  2.3 
days  to  45.5  days.  The  closest  orbit  (-22,0)  made  six  orbits  (2.3  days) 
before  being  perturbed  into  the  moon’s  surface.  A  trajectory  plot  for  this 
orbit  is  shown  in  Figure  9.  If  the  trajectory  is  continued,  ignoring  the  sur¬ 
face  impact,  several  mote  orbits  are  made  before  escaping  into  an  eccentric 
orbit  about  Mars.  Figure  10  shows  a  typical  example  of  an  extended  tra¬ 
jectory.  Allowing  the  trajectory  to  continue  after  a  known  collision  gives 
the  surface  of  section  plots  more  data  points  which  may  begin  to  show 
some  useful  pattern. 


Table  1 

Momenta  Calculations  for  Phobos 


Initial  Condition 
(-22,0) 
(-68,0) 
(-93,0) 
(-112,0) 
(-129,0) 
(-143,0) 
(-156,0) 
(-169,0) 
(-179,0) 
(-332,0) 
(-708,0) 


£1 

2.174818341961 
2.185520332825 
2.191359656599 
2.195808491525 
2.199797095646 
2.203087580257 
2.206147703997 
2.209212347677 
2.21 1572853458 
2.248029772709 
2.340508386779 
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Table  2 

Momenta  Calculations  for  Deimos 


Initial  Condition 


2x 


(-10,0) 

1.352473639530 

(-13,0) 

1.352646614364 

(-18,0) 

1.352934955024 

(-23,0) 

1.353223357308 

(-31,0) 

1.353684929236 

(-41,0) 

1.354262116378 

(-48,0) 

1.354666294454 

(-56,0) 

1.355128360725 

(-75,0) 

1.356226403753 

(-109,0) 

1.358193561438 

(-135,0) 

1.359699802981 

(-156,0) 

1.359699802981 

(-233,0) 

1.365392424480 

(-340,0) 

1.371635636656 

(-418,0) 

1.376205272893 
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All  five  inner  periodic  orbits  followed  this  same  typical  pattern  and 
collided  with  Phobos.  A  typical  surface  of  section  plot  is  shown  in  Figure 
1 1  for  the  (-22,0)  orbit  and,  ignoring  the  surface  impact,  shows  a  non-peri¬ 
odic  orbit  which  escapes  into  a  retrograde  orbit  about  Mars. 

Similarly,  the  six  outermost  orbits  remained  in  orbit  about  Phobos 
much  longer  than  the  closer  orbits  but  never  repeated  the  same  trajectory 
and  eventually  collided  with  the  moon  ranging  from  9  days  to  144  days. 
Again  the  surface  collision  was  ignored  in  order  to  see  if  any  pattern  in  the 
surface  of  section  plot  began  to  emerge.  A  typical  trajectory  is  shown  in 
Figure  12  for  the  (-143,0)  orbit  with  uneven  axes  for  visual  clarity.  Table 
3  summarizes  the  results  for  each  orbit. 

The  surface  of  section  plots  for  these  orbits  all  showed  the  same  pat¬ 
tern  of  non-periodic  orbits  which,  ignoring  the  surface  collision,  escape 
into  an  eccentric  orbit  about  Mars.  Figures  13  and  14  show  a  surface  of 
section  plot  for  the  (-143,0)  orbit  which  clearly  indicate  the  surface  impact 
(on  the  48th  day)  and  eventual  escape  into  an  eccentric,  prograde  orbit 
about  Mars.  The  same  pattern  is  seen  in  Figures  15  and  16  for  the  (-156,0) 
orbit  as  well  as  for  the  remaining  outermost  orbits. 

The  most  promising  case  was  the  (-332,0)  orbit  which  remained 
around  Phobos  the  longest  time  (144  days)  before  its  trajectory  collided 
with  the  surface.  Figures  17  and  18  clearly  show  the  collision  and  subse¬ 
quent  escape  of  the  spacecraft  to  a  prograde  orbit  about  Mars. 

With  the  exception  of  the  (-708,0)  orbit,  the  remaining  five  outer¬ 
most  periodic  orbits  remained  around  Phobos  well  over  45  days  before 
colliding  with  the  surface.  The  (-708,0)  orbit  is  the  furthest  periodic  orbit 
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Table  3 

Summary  of  Trajectory  Runs  for  Phobos 


Trajectory 

(-22,0) 

collision  at  2.3  days 

(-68,0) 

collision  at  10.7  days 

(-93,0) 

collision  at  26  days 

(-112,0) 

collision  at  28.6  days 

(-129,0) 

collision  at  45.5  days 

(-143,0) 

collision  at  48  days 

(-156,0) 

collision  at  50.3  days 

(-169,0) 

collision  at  69  days 

(-179,0) 

collision  at  71.3  days 

(-332,0) 

collision  at  144  days 

(-708,0) 

collision  at  9  days 

from  Phobos  found  by  Jansson  and  it  collided  with  the  surface  in  only  9 
days. 

With  only  the  eccentricity  of  Phobos  added  to  the  system  dynamics, 
aU  the  periodic  orbits  found  by  Jansson  became  unstable  and  collided  with 
the  moon.  Therefore,  no  other  perturbations  were  added  to  the  system  for 
study. 

Deimos 

Table  4  is  a  summary  of  the  results  for  the  orbits  investigated  around 
Deimos  with  only  the  eccentricity  added  to  the  system  dynamics.  Again, 
the  inner  orbits  are  of  most  interest. 

With  only  the  eccentricity  of  Deimos  added  to  the  system  dynamics, 
the  first  seven  inner  periodic  orbits  found  by  Jansson  either  collided  with 
the  surface  of  Deimos  or  escaped  into  an  eccentric  orbit  about  Mars. 

Figures  19  to  25  graphically  show  the  collisions  and  escapes  for  these  seven 
orbits  which  ranged  from  2.5  hours  to  9.9  days.  As  expected  for  the  inner 
orbits,  as  the  orbits  moved  away  from  Deimos  they  increased  their  orbital 
periods  and  improved  their  chances  of  remaining  in  orbit  longer.  Figure 
25  for  the  orbit  at  (-48,0)  clearly  shows  six  orbits,  with  a  closest  approach 
of  about  25  km,  before  escaping  into  a  retrograde  orbit  about  Mars.  The 
best  surface  of  section  plot  for  the  inner  orbits  is  shown  in  Figure  26  and 
displays  a  random  scattering  of  crossings  for  the  (-48,0)  orbit. 

The  eight  outermost  periodic  orbits  found  by  Jansson  all  behaved 
much  better  and  stayed  around  Deimos  for  the  entire  integration  period  of 
115  days  without  colliding  with  the  surface  or  escaping.  In  order  to  see  if 
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Table  4 

Summary  of  Trajectories  for  Deimos 


Initial  Condition  Trajecm 


(-10,0) 

Collision  at  2.5  hours 

(-13.0) 

Escaped  at  23.3  hours 

(-18.0) 

Escaped  at  2.1  days 

(-23.0) 

Escaped  at  2.3  days 

(-31.0) 

Collision  at  15  hours 

(-41.0) 

Collision  at  2.3  days 

(-48,0) 

Escaped  at  9.9  days 

(-56,0) 

Remained  for  578  days 

(-75,0) 

Remained  for  578  days 

(-109,0) 

Remained  for  578  days 

(-135,0) 

Remained  for  578  days 

(-156,0) 

Remained  for  578  days 

(-233,0) 

Remained  for  578  days 

(-340,0) 

Remained  for  578  days 

(-418,0) 

Remained  for  578  days 

any  pattern  developed  for  these  outer  orbits,  the  integration  period  was 
increased  to  578  days.  This  increased  the  integration  time  to  over  60 
minutes  for  each  run.  A  typical  trajectory  plot  for  the  seven  outermost 
orbits  is  shown  in  Figure  27  and  shows  an  oscillating  eccentric  orbit 
remaining  around  Deimos.  Based  on  the  appearance  of  the  trajectory  and 
the  length  of  the  integration  period  without  a  collision  or  escape,  these 
eight  outermost  orbits  are  considered  stable. 

Figure  28  shows  the  surface  of  section  plot  for  the  (-56,0)  orbit  and 
has  no  closed  circles  or  ellipses  which  would  indicate  that  there  are  no 
periodic  orbits.  However,  since  this  is  a  time  varying  Hamiltonian  system 
the  surface  of  section  technique  does  not  strictly  apply;  if  the  added  pertur¬ 
bation  is  small  enough,  one  could  expect  to  see  the  surface  of  section  points 
collected  within  some  region  or  boundary.  This  is  clearly  seen  in  Figure 
28. 

Figures  29  thru  37  show  more  examples  of  the  surface  of  section 
plots  for  these  seven  outermost  orbits  of  Deimos.  Although  not  as  random 
in  pattern  as  in  Figure  28,  they  do  appear  as  if  they  will  fill  in  the  same 
bounded  area  but  at  a  much  slower  pace.  Several  close-up  views  are  pre¬ 
sented  to  get  a  better  look  at  the  structure.  These  close-up  views  clearly 
show  there  are  no  closed  circles  or  ellipses  and  that  they  appear  as  if  they 
will  loop  back  and  continue  to  fill  in  the  same  bounded  area  as  in  Figure 
28. 

When  the  ^2  term  of  Mars’  geopotential  was  added  to  the  the  system 
dynamics  the  (-56,0)  orbit  also  went  unstable  and  collided  with  the  moon 
after  130  days.  The  last  seven  orbits  continued  to  remain  stable. 
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vn.  Conclusions 


All  integration  runs  were  made  with  just  the  eccentricity  of  the 
moons  added  to  the  system  dynamics.  Since  all  of  the  periodic  orbits  pre¬ 
viously  found  around  Phobos  became  unstable,  no  other  perturbations 
needed  to  be  introduced  to  the  system  to  prove  the  orbits  were  unrealistic. 
In  a  similar  manner,  the  seven  innermost  periodic  orbits  of  Deimos  which 
are  of  more  practical  use  also  became  unstable  with  just  the  eccentricity  of 
Deimos  added  to  the  system  dynamics. 

Phobos 

It  is  quite  clear  that  for  Phobos  all  the  periodic  orbits  found  by 
Jansson  were  based  on  the  assumption  that  the  orbit  of  Phobos  was  circular. 
However,  the  true  eccentricity  of  .015  proved  to  be  too  large  of  a  pertur¬ 
bation  for  the  assumption  to  hold  and  aU  the  orbits  became  unstable  and 
collided  with  the  surface  of  Phobos.  Even  though  the  orbits  are  not  peri¬ 
odic  in  nature,  any  spacecraft  attempting  to  rendezvous  with  another  planet 
will  have  the  capability  to  make  station  keeping  maneuvers.  With  the  range 
of  orbital  times  from  2.3  days  to  144  days,  there  are  plenty  of  opportuni¬ 
ties  for  these  type  of  corrections.  It  is  unfortunate  that  the  closest  two 
orbits  collide  with  Phobos  in  less  than  1 1  days  making  them  unsafe  as 
choices  for  a  paiking  orbit  even  with  station  keeping.  The  next  closest 
orbit  is  located  at  (-93,0)  and  gives  a  reasonable  safety  margin  to  make  a 
station-keeping  maneuver  before  colliding  with  the  surface  in  26  days. 

With  the  exception  of  the  (-708,0)  orbit,  the  safety  margins  increased  with 
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increasing  distance  from  Phobos,  allowing  more  flexibility  but  at  the  same 
time  reducing  its  desirability  as  a  parking  orbit. 

The  (-708,0)  orbit  is  the  furthest  away  from  Phobos  and,  as  seen  in 
Table  3,  did  not  follow  the  expected  trend  of  increasing  orbital  time  before 
a  collision  occurred.  This  is  may  be  due  to  the  satellite’s  attempt  to  enter  a 
true  orbit  around  Phobos  instead  of  an  eccentric  orbit  around  Mars,  that 
stays  close  to  Phobos.  This  attempt  is  doomed  by  the  fixed  initial  momenta 
value  based  on  Phobos’  revolution  rate  around  Mars. 

Deimos 

Because  of  the  much  smaller  eccentricity  of  Deimos  than  Phobos  and 
the  much  larger  distance  from  Mars,  the  orbits  around  Deimos  were 
expected  to  behave  better  than  those  of  Phobos.  The  seven  innermost 
orbits  still  did  not  handle  the  eccentricity  perturbation  very  well  and  either 
collided  or  escaped  rather  rapidly.  The  eighth  orbit  (-56,0)  went  unstable 
with  the  addition  of  the  geopotential  term  of  Mars  while  the  seven  outer¬ 
most  orbits  did  settle  down  into  an  acceptable  stable  condition. 

Of  the  seven  innermost  orbits,  four  escaped  from  Deimos  and  there¬ 
fore  could  be  used  safely  without  fear  of  imminent  collision  with  the  sur¬ 
face.  Of  course  station  keeping  maneuvers  could  easily  be  made  to  keep 
the  spacecraft  around  indefinitely.  The  three  collision  orbits  ranged  in 
orbital  time  from  2.5  hours  to  9.9  days  and  may  not  be  safe  choices  for 
parking  orbits  even  with  station  keeping. 

Although  the  seven  outermost  orbits  of  Deimos  did  appear  to  be 
stable  they  are  not  considered  very  useful  for  a  Mars  mission  due  to  their 
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much  larger  distance  from  the  Martian  surface  and  larger  distances  from 
Deimos  itself  which  increases  its  orbital  period  and  decreases  its  usefulness 
as  a  way-station.  AdditionaUy,  the  much  larger  distance  from  Deimos  to 
Mars  realistically  means  there  will  be  a  negligible  effect  from  the  J2  term 
of  the  geopotential,  and  the  perturbing  force  of  the  sun  can  realistically  be 
considered  negligible.  For  these  reasons,  and  the  fact  that  Phobos  is  con¬ 
sidered  a  more  useful  choice  for  a  way-station,  the  perturbing  force  of  the 
sun  was  not  added  to  the  dynamics  to  fully  test  the  seven  outermost  orbits 
of  Deimos. 

These  seven  outermost  orbits  did  remain  around  Deimos  for  the 
entire  integration  period  of  578  days,  oscillating  along  the  Y  axis  of 
Deimos,  as  shown  in  the  figures  described  earlier.  The  surface  of  section 
plots  indicate  a  boundary  constraining  the  points  and  was  expected  for 
small  perturbations.  Any  one  of  these  seven  orbits  could  be  used  as  a 
parking  orbit  for  Deimos  without  the  need  for  any  station  keeping  maneu¬ 
vers.  Realistically,  only  the  orbits  under  100  km  would  most  likely  be 
considered  as  serious  candidates  for  paiicing  orbits,  and  only  if  Deimos  is 
seriously  considered  as  a  way-station  over  Phobos. 


Appendix  A  Problem  Parameters  (4:67) 


Mqohs:  PhbbQS 

Axis  Lengths  (km)  x  =  13.4 

y  =  11.2 
z  =  9.2 


Moments  of  Inertia  (  kg*km2  ) 

Ixx  =  42.016  *  M  moon 
lyy  =  52.840  *  M  moon 
Izz  =  61.000  *  Mmoon 


Average  Semi-Major  Axis  Radius  (km) 

a  =  9,378 


Rotation  Rate  (rad/sec) 


O  =  0.000228 


Gravitational  Parameter  (km^/sec^) 

GMmoon  =s  0.00066 

Orbital  Eccentricity  (12:6.28) 

e  =  0.015 


Mars: 

Gravitational  Parameter  (km^/.sec^)  GMmi 

Equatorial  Radius  (km)  (12:6.3)  Rmbo 

Geopotential  Coefficient  (12:6.7)  h  = 


PgimbS 

x  =  7.5 
y  =  6.1 
z  =  5.2 


Ixx  =  12.850  *  M  moon 
lyy  =  16.658  *  M  moon 
Izz  =  18.692  *  Mnaoon 


a  =  23,459 


Q  =  0.0000576 


GMmoon  =  0.000088 


e  =  0.00052 


«  =  42828.32 
=  3397.2 
.0019557 
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Figure  10.  Extended  Trajectory  for  Orbit  (-22,0) 
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Figure  11.  Surface  of  Section  Plot  for  Orbit  (-22,0) 
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Figure  14.  Surface  of  Section  Plot  for  Orbit  (-143,0) 
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Figure  15.  Surface  of  Section  Plot  for  Orbit  (-156,0) 
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Figure  18.  Surface  of  Section  Plot  for  Orbit  (-332,0) 
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Figure  19.  Collision  Trajectory  for  Orbit  (-10,0) 
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Figure  20.  Escape  Trajectory  for  Orbit  (-13,0) 


67 


so  -30  -10  10  30 

X  (km) 


Figure  21.  Escape  Trajectory  for  Orbit  (-18, 
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Figure  22.  Escape  Trajectory  for  Orbit  (-23,0) 
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Figure  23.  Collision  Trajectory  for  Orbit  (-31,0) 
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Figure  24.  Collision  Trajectory  for  Orbit  (-41,0) 
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Figure  25.  Escape  Trajectory  for  Orbit  (-48,0) 
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Figure  26.  Surface  of  Section  Plot  for  Orbit  (-48,0) 
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Figure  27.  Trajectory  Plot  for  Orbit  (-56,0) 
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Figure  28.  Surface  of  Section  Plot  for  Orbit  (-56,0) 
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Figure  29.  Surface  of  Section  Plot  for  Orbit  (-75,0) 
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Figure  30.  Surface  of  Section  Plot  for  Orbit  (-75,0) 
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Figure  31.  Surface  of  Section  Plot  for  Orbit  (-109,0) 
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Figure  32.  Surface  of  Section  Plot  for  Orbit  (-109,0) 


79 


Figure  37.  Surface  of  Section  Plot  for  Orbit  (-233,0) 


84 


Bibliography 


1.  Croswell,  K.  “Target:  Phobos.”  Ad  Astra.  1:  17-21  (January  1989). 

2.  Lemer,  E.  J.  “Phobos:  Gains  from  a  Loss,”  Aerospace  America.28: 
42-46  (February  1990). 

3.  Goldman,  S.  “The  Legacy  of  Phobos  2,”  Skv  and  Telescope.  156-159 
(February  1990). 

4.  Jansson,  Capt  S.  W.  Stable  Orbits  about  the  Martian  Moons.  MS 
Thesis  AFnyGA/ENY/89D-3.  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wright-Patterson  AFB  OH, 

December  1989. 

5.  Brouwer,  D.  and  G.  Clemence.  Methods  of  Celestrial  Mechanics. 
New  York:  Academic  Press,  1961. 

6.  Meirovitch,  L.  Methods  of  Analytical  Dynamics.  New  York: 
McGraw-Hill  Book  Company,  1970. 

7.  Likens,  P.W,  Elements  of  Engineering  Mechanics.  New  York: 
McGraw-Hill  Book  Company,  1973. 

8.  Wiesel,  W.  E.  Adyanced  Astrodynamics.  Class  Lecture  Notes  for 
Mech  636,  School  of  Engineering,  Air  Force  Institute  of  Technology 
(AU),  Wright-Patterson  AFB,  OH,  December  1989. 

9.  Szebehely,  Victor.  Theory  of  Orbits.  New  York:  Academic  Press 
Inc.,  1967. 

10.  Cheney,  W.  and  D.  Kincaid.  Numerical  Mathematics  and  Computing. 
(Second  Edition).  Monterey,  CA.:  Brooks/Cole  Publishing 
Company,  1985. 

11.  Bates,  R.,  D.  Mueller  and  J.  White.  Fundamentals  of  Astrodynamics. 
New  Yoric:  Doyer  Publications,  Inc.,  1971. 


85 


12.  Smith,  R.  and  G.  West.  Space  and  Planetary  Environment  Criteria 
Guidelines  for  Use  in  Space  Vehicle  Developmem,  1982  Revision: 
Vol  1.  NASA  Technical  Memorandum  82478.  Marshall  Space 
Flight  Center,  Alabama,  1983. 


86 


Captain  Samuel  R.  Bryant  received  the  degree  of  Bachelor  of  Science 
in  Aerospace  Engineering  from  Mississippi  State  University  in  1984  and 
was  inducted  into  the  National  Honor  Society  for  Aerospace  Engineers.  In 
1985,  he  received  his  Air  Force  Commission  as  a  distinguished  graduate 
from  Officer  Training  School  in  San  Antonio,  Texas.  He  was  then  assigned 
to  HQ  Space  Division,  Los  Angeles  AFB,  CA.,  where  he  served  as  Mission 
Integration  Manager  for  the  Defense  Satellite  Communications  System 
(DSCS)  Phase  HI  satellite  program.  During  the  next  four  years.  Captain 
Bryant  was  responsible  for  coordinating  all  mission  integration  efforts  for 
transporting,  mating,  testing,  and  launching  the  DSCS  m  satellite  with  its 
assigned  launch  vehicle  which  included  the  Titan  34-D,  Titan  FV,  Space 
Shuttle,  and  the  newly  acquired  Atlas  11.  He  then  entered  the  Graduate 
School  of  Engineering  at  the  Air  Force  Institute  of  Technology  in  May 
1989. 

Captain  Bryant  is  a  member  of  Sigma  Gamma  Tau,  AlAA,  and  the 

AFA. 


Permanent  Address: 


P.O.  Box  55 

Port  Bane,  LA.  70577 


I  REPORT  DOCUMENTATION  PAGE  ( 

1.  AGENCY  USE  ONEV  2  RETORT  O.ME  3.  RE‘'OKr  TYPE  AND  DATES  COVERED 

December  90  Final 

4.  TITLE  AND  SUBTITLE 

Realistic  Orbits  About  the  Martian  Moons 

5.  FUNDING  NUMBERS 

6.  AUTHOR(S) 

Samuel  R.  Bryant,  Captain,  USAF 

7.  PERrORMiNG  OKGANIZATIC.N  NAME(S;  AND  ADDRESSIES) 

School  Of  Engineering 

Air  Force  Institute  of  Technology 

£  PERFORMING  ORGANIZATION 

REPORT  NUMBER  j 

AFIT/GA/ENY/90D-1  | 

1 

9.  SPONSORING  MONITORING  AGENCY  NAME'S)  AND  ADDR£SS(ES) 

10.  SPONSORING  MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

12d.  distribution  code 

13.  abstract  (Maximurri  ^’00  if-'O'ClS'- 

Orbits  about  the  Martian  moons,  Phobos  and  Deimos,  previously  found  to  be  stable  by 
the  Poincare  surface  of  section  technique  in  a  restricted  three-body  system,  were 
investigated  with  non-zero  eccentricity  of  the  moons  added  to  the  dynamics  of  the 
planar  four-body  system  (Mars,  moon,  sun,  artificial  satellite).  Although  not 
strictly  applicable  to  this  case,  the  surface  of  section  technique  was  used  to 
verify  whether  or  not  the  previously  found  stable  orbits  were  realistic  when  a  real 
world  perturbation  was  added  to  the  system  dynamics.  The  surface  of  section  technique 
involves  the  numerical  integration  of  several  orbits  with  the  same  initial  value  of 
the  Hamiltonian  based  on  matching  the  angular  rates  of  the  orbit  and  moon  with  respect 
to  Mars.  Apoapsis  and  periapsis  points  of  the  orbits  were  plotted  in  the  two- 
dimensional  configuration  space.  Stable  orbits  were  found  when  the  points  formed 
regions  which  are  bounded;  chaotic  orbits  were  indicated  by  random  points. 


14.  SUBJECT  TERMS 

Mars  Phobos  Deimos  Poincare  Surface 

of  Section 

IS  NUMBER  OF  PAGES 

99 

Orbits 

16  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

16.  security  Classification 

OF  THIS  PAGE 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

20  limitation  Of  abstract 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 

UNLIMITED 

NS^.  ’5iG.u'.:80-5500  .'98  .  89 


